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The existence of a quantum percolation threshold p q < 1 in the 2D quantum site-percolation 
problem has been a controversial issue for a long time. By means of a highly efficient Chebyshev 
expansion technique we investigate numerically the time evolution of particle states on finite disor- 
dered square lattices with system sizes not reachable up to now. After a careful finite-size scaling, 
our results for the particle's recurrence probability and the distribution function of the local particle 
density give evidence that indeed extended states exist in the 2D percolation model for p < 1. 
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I. INTRODUCTION 

For classical percolation above a critical concentration 
of accessible sites, the so-called percolation threshold p Cl 
a spanning cluster exists, allowing for transport. In the 
quantum case, however, the existence of the spanning 
cluster does not guarantee an extended wave function. 
As for the Anderson 1 and binary alloy 2 disorder models 
scattering and interference effects at the irregular bound- 
aries of the cluster may lead to localization of the quan- 
tum particle, i.e. absence of diffusion. Increasing the con- 
centration of accessible sites p above p Cl the first occur- 
rence of extended states defines a quantum percolation 
threshold, p q . For the Anderson model, below a critical 
disorder strength, a single mobility edge separates bands 
of localized and extended states. By contrast, for the bi- 
nary alloy and three-dimensional (3D) percolation model 
a sequence of "mobility edges" exists 3 . It is clear that 
Pc < Pq < 1 and the fundamental question is, of course, 
whether p q equals one of these boundaries. 

agree on 
3jD ranging from 0.4 
to 0.5 for site percolation on a simple cubic lattice. In 2D, 
the situation is less clear. Here the physical community 
is evenly divided: one group^ 7 -^ favors p 2D — 1 while 
another group^^^i 3 -^^^ claims that p 2 q D < 1. The 
most striking argument against p 2D < 1 comes from one- 
parameter scaling theory 17 , according to which arbitrary 
small disorder always leads to localization in 2D. Never- 
theless, there are hints for p 2D < 1. In this regard band 
center states seem to be of particular importanc e 18 ^ 9 . 
Whether E = states on a 2D bipartite depleted square 
lattice are a possible "let out clause"—^ for the one- 
parameter scaling theory result is an open question. 

Renewed interest in 2D quantum percolation came up 
in connection with the unusual transport properties of 
novel materials. For example, the metal-insulator transi- 
tion in perovskite manganite films and the related colos- 
sal magnetoresistance effect seems to be inherently per- 
colativ o 22 i 23 . Quantum percolation might be of impor- 
tance for the experimentally observed metallic behav- 
ior of dilute weakly disordered Si-MOSFETs 24 . More- 
over, quite recently, a quantum percolation scenario 



For the 3D case, results in the literature 3 -^ 
Pc D < Pq D < 1? with estimates for p\ 



has been proposed for minimal conductivity in undoped 
graphene 25 . All this gives additional motivation for the 
rather fundamental study of 2D quantum percolation 
presented in this work. 



II. MODEL AND METHODS 

We start from a tight-binding Hamiltonian of non- 
interacting spinless fermions in Wannier representation 
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with uniform hopping t between nearest neighbors (ij) 
on a finite square lattice with N = L 2 sites and periodic 
boundary conditions. The on-site energies are subject 
to the bimodal distribution 



p[ei] =pS(ci - ca) + (1 -p)S(ei - e B ) • 



(2) 



The quantum site-percolation Hamiltonian then is ob- 
tained in the limit |e# — ea\ — > oo. Without loss of gen- 
erality we choose ca = 0. In this situation the fermions 
move on a random assembly of A-lattice points. Depend- 
ing on the largest cluster of TV connected A-sites may 
span the entire lattice. Note that the dynamics of our fi- 
nite system is governed by two different time scales. The 
time scale of a hopping process is given by the inverse 
of the hopping matrix element, To = 1/t. To account for 
the system size, we define a characteristic time To = iVYo, 
which, in principle, is necessary to visit each site of the 
spanning cluster once. 

In order to address the problem of localization in quan- 
tum percolation, we determine, for p > p c = 0.593, the 
recurrence probability of a particle to a given site, Pr(t), 
which for r — » oo may serve as a criterion for Anderson 
localizatio n 26 ! 27 . While for extended states on the span- 
ning cluster Pr(t — > oo) = 1/N scales to zero in the 
thermodynamic limit, for localized states Pr(t) remains 
finite as TV — > oo. 

Concept ionally, in a first step, we track the time evo- 
lution of an initially localized state (which, of course, in 
general is not an eigenstate) by calculating the modulus 
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square of the wave function at the particle's starting po- 
sition as a function of time r. To this end we expand 
the time evolution operator U(r) = e~ lHr in Cheby- 
shev polynomials 2 ^ 29 ! 30 , where the Hamiltonian has to 
be rescaled to the definition interval of the Chebyshev 
polynomials T k: [—1,1]. With H = aH + 6 we obtain 
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The expansion coefficients c k are given by 
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where J k denotes the Bessel function of the first kind. To 
calculate the evolution of a state |VK r )) fr° m one time 
grid point to the next, \ip(r + At)) = J7(Ar)|^(r)), we 
have to accumulate the c^-weighted contributions of the 
different orders \v k ) — T k (H)\ip(r)). As the coefficients 
c k (a At) depend on the time step but not on time explic- 
itly, we need to calculate them only once. The contri- 
butions \v k ) can then be calculated iteratively using the 
recurrence relation of the Chebyshev polynomials 



\vk+i) = 2H\v k ) - \v k -i) 



(5) 



where \v\) = H\vo) and \vq) = \iJj(t)). Due to the fast 
asymptotic decay of the Bessel functions 
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the coefficients c k vanish rapidly above a certain expan- 
sion order. As a result from the infinite series only a fi- 
nite number of M terms needs to be taken into account. 
Clearly M depends on the time step used. Truncating the 
series for large ar > 10 3 at M ~ 1.5ar, the \c k \ < 10 -16 
for all discarded terms. For even larger values of ar, 
the necessary value of M/(ar) is even reduced and ap- 
proaches M ~ ar for r — > oo. Thus, as compared to 
the standard Crank-Nicolson algorithm, the Chebyshev 
expansion permits the use of a considerably larger time 
step to achieve the same accuracy, i.e. allows for a very 
efficient calculation of the time evolution of a given state. 

In a second step, we employ the so-called local distri- 
bution approac h 31 ! 32 . In the theoretical investigation of 
disordered systems it turns out that the distribution of 
random quantities is central. While all characteristics 
of a certain material are determined by the distribution 
p[ci], each actual sample constitutes only one particular 
realization {q}. At each randomly chosen site i of such 
a particular sample, we observe different local environ- 
ments. That is disorder breaks translational invar iance 
and we have to focus on site-dependent quantities like 
the local density of states (LDOS) 33 at site i 
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Given an energy E, the LDOS can be efficiently calcu- 
lated by means of the kernel polynomial method 34 . Con- 
sidering the LDOS, a well established criterion for local- 
ization is the following. Probing different sites in the 
sample and recording the values of pi gives the proba- 
bility distribution f[pi\. To alleviate the problem of sta- 
tistical noise one usually considers instead of f[pi] the 
distribution function 
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In the thermodynamic limit both f[pi] and F[pi] are 
independent of the actual realization {e^} (i.e. self- 
averaging), but solely depend on p[ci\. Therefore F[pi] 
characterizes H(p[ci\), not only H({ei}). That means, at 
the level of distributions, we have restored translational 
invariance. For an extended state the amplitude of the 
wave function is more or less uniform and f[pi] is sharply 
peaked and symmetric around its mean value {pi). For 
localized states pi strongly fluctuates throughout the lat- 
tice, giving rise to a very asymmetric f[pi] with a long 
tail and {pi) — > 0. Consequently the distribution func- 
tion F[pi] steeply rises for extended states whereas for 
localized states the increase extends over several orders 
of magnitude. These arguments also hold for the (time- 
dependent) particle density at site i 
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if we take additional care of some aspects. Evolving an 
arbitrary initial state \ip(0)) in time, we only have access 
to 7ii(r), containing contributions of the whole spectrum. 
Calculating the time evolution of an initial state we have 
access to rii on the whole cluster, which gives a much 
better statistics. On the other hand, the chosen initial 
conditions introduce an additional dependency. Espe- 
cially for localized states the local environment of the 
starting position will influence F[nJ. So for n^, in ad- 
dition to taking the thermodynamic limit, we have to 
examine different initial conditions. While the LDOS 
approach allows for an energy resolved examination of 
the localization problem, rii gives no information about 
p q (E) because it contains contributions from all E. But, 
starting from p c and increasing p, we detect at a certain 
P — Pq the first occurrence of extended states somewhere 
in the spectrum. 



III. RESULTS 

Let us now discuss the outcome of our numerical in- 
vestigations. Figure [1] compares the time evolution of an 
initially localized state on the spanning cluster for low 
(left panels) and high (right panels) concentration of A 
sites, for which qualitatively different behaviors emerge. 
For small p = 0.65 > p c , the initial state spreads within 
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FIG. 1: (Color online) Time evolution of the normalized local 
particle density Nm(r) for an initially localized state on the 
spanning cluster of a N — 512 2 square lattice for p = 0.65 
(left) and p = 0.90 (right). Due to normalization, for an 
extended state evenly spread over all sites of the spanning 
cluster this quantity is equal to unity. 



a short time over a finite region of the spanning cluster. 
But also for very long times, r > 100Tb, this extension 
does not change significantly anymore. For larger p, the 
spreading is even faster. In contrast to the previous case, 
for p = 0.90 the state is not confined to some finite re- 
gion, but rii is transfered to the whole cluster. Since 
the initial state is a superposition of eigenstates from the 
whole spectrum, it also contains contributions from local- 
ized states, for instance near the band edges. Those are 
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FIG. 2: (Color online) Time evolution of the distribution func- 
tion of the local particle density, F[Nm(r)]. The curves corre- 
spond to the data in Fig. [T] supplemented by the distribution 
functions at r = (dot-double-dashed curve). 

reflected in the darker spots in the vicinity of the initial 
position, which persist there, even for very long times. 
This signature of localized states we already know from 
the case of small p. The essential new feature is that for 
large p some eigenstates exist, which are not localized, 
i.e. spread over the whole cluster. 

To go beyond this simple visualization and account for 
a quantitative description, in Fig. [2] we show the cor- 
responding local particle density distribution functions 
F[rii]. The different behavior of localized and extended 
states becomes obvious as the time increases. Since these 
results were obtained for a particular sample and a cer- 
tain starting position it is necessary to ask how strong 
these findings depend on the initial conditions, and if 
they are representative for the underlying p[ei\. 

To answer this question we examine systematically the 
influence of the initial conditions on the behavior of the 
state. In this respect the geometry and energy of the 
initial state, the starting position on the cluster as well 
as the particular cluster realization {e^} should be of im- 
portance. A fingerprint in which all these aspects come 
into play is the LDOS of the initial state. 

In Fig.[3]we examine the influence of these issues one by 
one. Fixing a particular cluster realization and starting 
position of the wave function, we may construct several 
initial states with the same total energy. This can be 
done by choosing a different number of sites which ini- 
tially have non- vanishing amplitudes. For a state with 
finite amplitude only on two neighboring sites, a and 
\f \ — a 2 , respectively, we get E = 2ta\]\ — a 2 and we 
may continuously tune E G [— t, t] by choice of a (cf. 
Fig. [3] a and c). These energies may also be constructed 
for more complicated initial states, as e.g. for the one 
shown in Fig. [3]b, with non- vanishing amplitudes on six 
adjacent sites. The more complicated these structures 
get and the more sites are involved, the larger energies 
of the initial state are possible, e.g. for the configuration 
of Fig. [3]b up to (1 + V2)t. As stressed above, these en- 
ergies are not eigenenergies of H, i.e. the initial state is 
a superposition of eigenstates from the whole spectrum. 
A quantitative characterization of the amount of differ- 



4 



ent contributions can be obtained by the LDOS of the 
initial state. Changing the starting position (Fig. Eld) or 
the cluster realization (Fig. [3]e) has a similar effect. Es- 
sentially, we see the same behavior as for the changes 
in Figs. [3] a-c. Most notably, although the LDOS at 
r = and the evolution of the states differs in detail, 
on a coarse grained scale they behave similarly. This 
is also corroborated by the behavior of the distribution 
function F[n^(r)] for sufficient long times r. Despite mi- 
nor differences between the five curves, they agree on 
showing the characteristics of a localized state (Fig. [3] f, 
solid lines). The more a state is localized, i.e. the lower 
the more the local structure of the spanning cluster in- 
fluences its evolution and thus F[m(r)]. At large p, the 
dependency of an extended state on the above criteria 
is much weaker, as in this case the local environment of 
the starting position is less important (Fig. [3] f, dashed 
lines). Overall, the characteristics of the time evolution 
is mainly determined by the cluster structure, the initial 
state has only minor influence. This holds as the random 
nature of the cluster guarantees a similar structure above 
a certain scale. Finally one may ask if the chosen bound- 
ary conditions or the linear extension (odd/even) of the 
lattice influences the dynamics. The latter distinction 
has been shown to be important for a bond percolation 
model because of the bipartiteness of the underlying lat- 
tice^. In our case, however, we did not see an influence 
on the results. 

In view of a finite-size scaling of the numerical data we 
first have to ensure that the obtained distribution func- 
tion has already become quasi-stationary. To check this, 
we calculate the L2 norm of the difference between two 
distribution functions at different times. Due to fluctu- 
ations, we cannot expect this quantity to vanish com- 
pletely even for large times. For r > r* « 0.1Tb, how- 
ever, we have reached the quasi-stationary regime for all 
considered system sizes (cf. Fig. 2]). At this time, the 
wave function has reached its maximum extension and 
the further development is governed by amplitude fluc- 
tuations from site to site. In this regime, the difference 
between the distribution functions does not depend on 
time anymore (inset of Fig.Hj), which is a clear indication 
of random fluctuations without additional drift motion. 
This allows for a determination of the quasi-stationary 
distribution function together with an error estimation 
and enables us to compare different system sizes. In 
view of computation time, the scaling r* ~ To ~ N 
together with the linear dependency on the dimension of 
the Hilbert space for the Chebyshev expansion leaves us 
with a scaling ~ (pN) 2 ~ p 2 L A . 

As soon as we have quasi-stationary distribution func- 
tions, their dependency on the system size may be ex- 
ploited for a quantitative distinction between localized 
and extended states (Fig. [5]). For the two different occu- 
pation probabilities, we observe completely different be- 
havior: while for p = 0.65 the distribution function shifts 
towards smaller values on increasing the system size, for 
p = 0.90 it is not affected by the change of system size at 
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FIG. 3: (Color online) Influence of the initial conditions on 
the evolution of the wave function on a N = 512 2 system 
at p = 0.65. The colormap is the same as the one in Fig(T] 
Panels a-d show the same cluster realization, with identical 
starting positions for a-c. All states have E — 0.5t, except 
for the case shown in panel c, where E — t. The insets give a 
magnification of rii(0). The corresponding LDOS of the initial 
states are given below each panel. In panel f the solid lines 
show the distribution function F[rii(r)] for r = 15.4Tb. For 
comparison the dashed lines indicate F[m(r)] for p = 0.90 
subject to the same differences in the initial conditions. 
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FIG. 4: (Color online) L2 norm of the difference between 
the distribution function F[rii(r)] at time r and the quasi- 
stationary distribution F[ni(r*)] for different system sizes. 
The inset shows the dependency of this quantity on the time 
difference for the N = 4096 2 system in the quasi-stationary 
regime. Note the different axis scales. 
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FIG. 5: (Color online) Finite-size scaling of F[Nni(r*)] in the 
quasi-stationary regime. 



tide states for the 2D quantum site-percolation model. 
The local particle density contains contributions of eigen- 
states from the whole energy spectrum and is therefore 
an ideal tool to globally investigate the localization prop- 
erties of the system. Examining the corresponding dis- 
tribution function allows for a distinction between local- 
ized and extended states. Above a certain concentration 
of accessible sites the shape of the distribution function 
becomes independent of the system size. This gives evi- 
dence that there are (some) extended states in the spec- 
trum. Having access to much larger systems than any 
others previously in the literature, we conclude from our 
data, supplemented by a careful finite-size scaling, that 



„2D < L 



all. The latter is clearly the behavior one would expect 
for an extended state. 
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